module bimocq_common_mod

  implicit none

contains

real function limit_coord(xmin, xmax, xrange, x) result(res)

  real, intent(in) :: xmin
  real, intent(in) :: xmax
  real, intent(in) :: xrange
  real, intent(in) :: x

  res = x
  if (res < xmin) then
    do while (res < xmin)
      res = res + xrange
    end do
  else if (res > xmax) then
    do while (res > xmax)
      res = res - xrange
    end do
  end if

  end function limit_coord

  subroutine interp1d(method, ns, nx, x, f, x0, f0)

    integer, intent(in) :: method
    integer, intent(in) :: ns
    integer, intent(in) :: nx
    real, intent(in) :: x(1-ns:nx+ns)
    real, intent(in) :: f(1-ns:nx+ns)
    real, intent(in) :: x0
    real, intent(out) :: f0

    integer n, i(4), p, q
    real dx, wx(4), x1, x2

    dx = x(2) - x(1)
    select case (method)
    case (1)
      n = 2
    case (2)
      n = 3
    case (3)
      n = 4
    end select
    ! Set indices of involved grids.
    i(1) = int(x0 / dx) + 1 - n / 2 + 1
    do p = 1, n - 1
      i(p+1) = i(p) + 1
    end do
    ! Calculate weights.
    do p = 1, n
      x1 = x(i(p))
      wx(p) = 1.0
      do q = 1, n
        if (p == q) cycle
        x2 = x(i(q))
        wx(p) = wx(p) * (x0 - x2) / (x1 - x2)
      end do
    end do
    ! Evaluate formula.
    f0 = 0.0
    do p = 1, n
      f0 = f0 + wx(p) * f(i(p))
    end do

  end subroutine interp1d

  subroutine interp2d(method, ns, nx, ny, x, y, f, x0, y0, f0)

    integer, intent(in) :: method
    integer, intent(in) :: ns
    integer, intent(in) :: nx
    integer, intent(in) :: ny
    real, intent(in) :: x(1-ns:nx+ns)
    real, intent(in) :: y(1-ns:ny+ns)
    real, intent(in) :: f(1-ns:nx+ns,1-ns:ny+ns)
    real, intent(in) :: x0, y0
    real, intent(out) :: f0

    integer n, i(4), j(4), p, q
    real dx, dy, wx(4), wy(4), x1, y1, x2, y2

    dx = x(2) - x(1)
    dy = y(2) - y(1)
    select case (method)
    case (1)
      n = 2
    case (2)
      n = 3
    case (3)
      n = 4
    end select
    ! Set indices of involved grids.
    i(1) = int(x0 / dx) + 1 - n / 2 + 1
    j(1) = int(y0 / dy) + 1 - n / 2 + 1
    do p = 1, n - 1
      i(p+1) = i(p) + 1
      j(p+1) = j(p) + 1
    end do
    ! Calculate weights.
    do p = 1, n
      x1 = x(i(p))
      y1 = y(j(p))
      wx(p) = 1.0
      wy(p) = 1.0
      do q = 1, n
        if (p == q) cycle
        x2 = x(i(q))
        y2 = y(j(q))
        wx(p) = wx(p) * (x0 - x2) / (x1 - x2)
        wy(p) = wy(p) * (y0 - y2) / (y1 - y2)
      end do
    end do
    ! Evaluate formula.
    f0 = 0.0
    do p = 1, n
      do q = 1, n
        f0 = f0 + wx(p) * wy(q) * f(i(p), j(q))
      end do
    end do

  end subroutine interp2d

end module bimocq_common_mod